##########################################################################################

library(dplyr)
library(data.table)
library(optparse)
library(ggplot2)
library(Seurat)
library(monocle)
library(ggpubr)

##########################################################################################
option_list <- list(
    make_option(c("--input_im_file"), type = "character"),
    make_option(c("--input_igc_file"), type = "character"),
    make_option(c("--input_dgc_file"), type = "character"),
    make_option(c("--igc_mono_file"), type = "character"),
    make_option(c("--dgc_mono_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--type"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 韩国的单细胞表达文件
    input_im_file <- "~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_TPM_ourData/MUC6_1.5/Scissor_STAD_MUC6_mutation.IM.CellRate.all.RData"
    input_igc_file <- "~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_TPM_ourData/MUC6_1.5/Scissor_STAD_MUC6_mutation.IGC.CellRate.all.RData"
    input_dgc_file <- "~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_TPM_ourData/MUC6_1.5/Scissor_STAD_MUC6_mutation.DGC.CellRate.all.RData"

    igc_mono_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/singleCell/njmu/IGC_DEG_MONO800.Rdata"
    dgc_mono_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/singleCell/njmu/DGC_DEG_MONO800.Rdata"

    gene <- "MUC6"
    type <- "all"
    ## 输出
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_TPM_ourData/MUC6_monocle"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_im_file <- opt$input_im_file
input_igc_file <- opt$input_igc_file
input_dgc_file <- opt$input_dgc_file
igc_mono_file <- opt$igc_mono_file
dgc_mono_file <- opt$dgc_mono_file
out_path <- opt$out_path
gene <- opt$gene
type <- opt$type

dir.create(out_path , recursive = T)

###########################################################################################

dat_im <- load(input_im_file)
dat_im <- sc_dataset

dat_igc <- load(input_igc_file)
dat_igc <- sc_dataset

dat_dgc <- load(input_dgc_file)
dat_dgc <- sc_dataset

dat_igc_time <- load(igc_mono_file)
dat_dgc_time <- load(dgc_mono_file)

###########################################################################################
scissor_1_cell <- c(
    names(which(dat_im$scissor==1)) ,
    names(which(dat_igc$scissor==1)) ,
    names(which(dat_dgc$scissor==1))
    )
scissor_2_cell <- c(
    names(which(dat_im$scissor==2)) ,
    names(which(dat_igc$scissor==2)) ,
    names(which(dat_dgc$scissor==2))
    )

IGC_DEG_MONO800_new <- IGC_DEG_MONO800
IGC_DEG_MONO800_new$scissor <- "Background"
IGC_DEG_MONO800_new$scissor[colnames(IGC_DEG_MONO800_new) %in% scissor_1_cell]<-"Scissor+"
IGC_DEG_MONO800_new$scissor[colnames(IGC_DEG_MONO800_new) %in% scissor_2_cell]<-"Scissor-"

DGC_DEG_MONO800_new <- DGC_DEG_MONO800
DGC_DEG_MONO800_new$scissor <- "Background"
DGC_DEG_MONO800_new$scissor[colnames(DGC_DEG_MONO800_new) %in% scissor_1_cell]<-"Scissor+"
DGC_DEG_MONO800_new$scissor[colnames(DGC_DEG_MONO800_new) %in% scissor_2_cell]<-"Scissor-"

col <- c('grey','indianred1','royalblue')
names(col) <- c("Background" , "Scissor+" , "Scissor-" )

###########################################################################################
## 描绘高变基因推测的时序在所有上皮细胞的情况
plotTrajectory <- function(ALL_cds_all = ALL_cds_all , out_path = out_path , type = type){

    ## 拟时序的时间分布
    image_name <- paste0( out_path , "/CellSubtype_trajectory.",type,".pdf" )
    plot2 <- plot_cell_trajectory(ALL_cds_all, color_by="celltype", cell_size=1)
    ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)

    image_name <- paste0( out_path , "/class_trajectory.",type,".pdf" )
    plot2 <- plot_cell_trajectory(ALL_cds_all, color_by="sample", cell_size=1)
    ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)

    image_name <- paste0( out_path , "/pseudotime_trajectory.",type,".pdf" )
    plot2 <- plot_cell_trajectory(ALL_cds_all, color_by="Pseudotime", cell_size=1)
    ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)

    image_name <- paste0( out_path , "/Patient_trajectory.",type,".pdf" )
    plot2 <- plot_cell_trajectory(ALL_cds_all, color_by="patient", cell_size=1)
    ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)

    image_name <- paste0( out_path , "/scissor_trajectory.",type,".pdf" )
    plot2 <- plot_cell_trajectory(ALL_cds_all, color_by="scissor", cell_size=1) + scale_color_manual(values=c(col))
    ggsave(file = image_name , plot = plot2,width = 6.5,height = 6)
}

type <- "IGC"
plotTrajectory( ALL_cds_all = IGC_DEG_MONO800_new , out_path = out_path , type = type)

type <- "DGC"
plotTrajectory( ALL_cds_all = DGC_DEG_MONO800_new , out_path = out_path , type = type)

###########################################################################################

tmp_data <- data.frame(table(dat_im@meta.data$scissor , dat_im@meta.data$celltype))
colnames(tmp_data) <- c("Scissor_Type" , "Cell_Type" , "Cells")
tmp_data <- tmp_data %>%
group_by( Cell_Type ) %>%
summarize( Scissor_Type = Scissor_Type , Cells = Cells , Cells_Rate = Cells/sum(Cells) )
cell_order <- names(table(dat_im$celltype)[order(table(dat_im$celltype) , decreasing=T)])

tmp_data$Scissor_Type <- as.character(tmp_data$Scissor_Type)
tmp_data$Scissor_Type <- ifelse( tmp_data$Scissor_Type==0 , "Background" , tmp_data$Scissor_Type )
tmp_data$Scissor_Type <- ifelse( tmp_data$Scissor_Type==1 , "Scissor+" , tmp_data$Scissor_Type )
tmp_data$Scissor_Type <- ifelse( tmp_data$Scissor_Type==2 , "Scissor-" , tmp_data$Scissor_Type )
tmp_data$Cell_Type <- factor( tmp_data$Cell_Type , levels = cell_order , order = T )

## 去除主细胞,该细胞数量太少且有MUC6突变定义不明确，以访误导
#tmp_data <- subset( tmp_data , Cell_Type != "Chief" )
col <- c('grey','indianred1','royalblue')
names(col) <- c("Background" , "Scissor+" , "Scissor-" )

tmp_data <- subset(tmp_data , Scissor_Type!="Background")
p1 <- ggplot(tmp_data,aes(x=Cell_Type,y=Cells,fill=factor(Scissor_Type))) +
    geom_bar(stat="identity") +
    ylab("Cell Counts") +
    xlab(NULL) +
    theme_bw() +
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='right',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        legend.text = element_text(size = 8,color="black",face='bold'),
        axis.text.x = element_text(size = 10,color="black",face='bold' , angle = 45, hjust = 1),
        axis.text.y = element_text(size = 8,color="black",face='bold'),
        axis.title.x = element_text(size = 10,color="black",face='bold'),
        axis.title.y = element_text(size = 10,color="black",face='bold'),
        strip.text.x = element_text(size = 15,color="black",face='bold'),
        axis.line = element_line(size = 0.5))  +
    scale_fill_manual(values=c(col)[2:3])

out_name <- paste0( out_path , "/Scissor_STAD_" , gene , "_mutation.IM.CellRate.all.pdf"  )
ggsave(file=out_name,plot=p1,width=6,height=6)